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We consider percolation and jamming transitions for particulate systems exposed to compression. For the 
systems built of particles interacting by purely repulsive forces in addition to friction and viscous damping, it 
is found that these transitions are influenced by a number of effects, and in particular by the compression rate. 
In a quasi-static limit, we find that for the considered type of interaction between the particles, percolation and 
jamming transitions coincide. For cohesive systems, however, or for any system exposed to even slow dynamics, 
the differences between the considered transitions are found and quantified. 
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I. INTRODUCTION 

The dense systems of particles interacting by either purely 
repulsive potentials, such as dry granular particles, or by both 
repulsive and attractive ones, such as wet granulates, appear 
virtually everywhere, from nature to a variety of applications 
bridging the scales from nano to macro. The structure of 
the force field by which the particles interact may be very 
complex, in particular on meso-scales where this force field 
is nonuniform and forms force networks. These networks 
are of relevance not only to granular systems, but to many 
other ones, such as foams and colloids. Their properties have 
been recently explored using a variety of different approaches, 
ranging from theoretical and computational ones based on ex¬ 
ploring local structure of force networks jj]|, networks type of 
approaches US , and topological methods di]. 

While percolation has been considered for dense particu¬ 
late systems CM, much more is known about static and or¬ 
dered lattice-based systems 030, for which two types of 
percolation are discussed - rigidity and connectivity percola¬ 
tion 130. However, lattice models do not account for non¬ 
linear effects at particle contacts, such as friction and viscous 
damping, or for dynamics, so it is unclear whether the results 
obtained for lattice systems apply to particulate ones 10 . For 
the latter, the connection between percolation (connectivity) 
and jamming (rigidity) transitions was discussed recently for 
both non-cohesive and cohesive frictionless systems, and it 
was found (for the systems considered) that these two tran¬ 
sitions in general differ JHI3]. However, these conclusions 
were reached by considering rather specific interaction mod¬ 
els (over-damped dynamics), and the question whether they 
hold in general, and whether they also follow from the models 
commonly used to simulate physical granular particles, is still 
open. 

In this paper, we discuss the relation between percolation 
and jamming for frictional and frictionless particles in two 
spatial dimensions, both with and without cohesion. We con¬ 
sider slowly compressed systems that go through percolation 
and jamming and discuss how these transitions depend on 
the system properties. The motivation for considering com¬ 
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pression is that it is a simple protocol that avoids the com¬ 
plexities associated with shear, and allow us to focus the dis¬ 
cussion. However, consideration of any dynamics, including 
compression, naturally leads to the questions related to the 
rate-dependence of the results, and, as we will see, to new 
insight into percolation and jamming transitions for evolving 
particulate systems. 

The paper is organized as follows. In Sec. [TT] we present 
the simulation techniques. In Sec.|III]we present out findings, 
first for purely repulsive systems in Sec. IHI A1 and then for 
cohesive ones in Sec. lIIIBl Section|Iy]is devoted to summary, 
conclusions, and future outlook. 

II. SIMULATIONS 

We perform discrete element simulations using a set of cir¬ 
cular particles confined in a square domain, using a slow- 
compression protocol 0,S, augmented by relaxation as de¬ 
scribed below. Initially, the system particles are placed on a 
square lattice and are given random velocities; we have ver¬ 
ified that the results are independent of the distribution and 
magnitude of these initial velocities. The discussion related 
to possible development of spatial order as the system is com¬ 
pressed can be found in @1, and the issue of spatial isotropy 
of the considered systems is considered later in the text. 

In our simulations gravity is not considered, and the diame¬ 
ters of the particles are chosen from a flat distribution of width 
r p . System particles are soft inelastic disks and interact via 
normal and tangential forces, including static friction, ,u (as 
in [3Bl). The particle-particle (and particle-wall) interactions 
include normal and tangential components. The normal force 
between particles i and j is 

F lj = k„xn - y n mx" j (1) 

r ij = Kj\, r jj = Tf-Tj, n - nj/nj 

where v" ■ is the relative normal velocity. The amount of com¬ 
pression is x = djj — rij, where dij = (dj + dj)/ 2, dj and 
dj are the diameters of the particles i and j. All quanti¬ 
ties are expressed using the average particle diameter, d aV e, 
as t he lengthsca le, the binary particle collision time x c = 
2 It \J d U ve /(2 gkf) as the time scale, and the average particle 
mass, m, as the mass scale, m is the reduced mass, k n (in 
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units of mg/dnye) is set to a value corresponding to photoe¬ 
lastic disks 111511 . and y„ is the damping coefficient Ilq]. The 
parameters entering the linear force model can be connected 
to physical properties (Young modulus, Poisson ratio) as de¬ 
scribed e.g. in 111611 . 

We implement the commonly used Cundall-Strack model 
for static friction JB, where a tangential spring is intro¬ 
duced between particles for each new contact that forms at 
time t = to. Due to the relative motion of the particles, the 
spring length, £, evolves as ^ = ff v' ; ( t') dt', where v' ■ = 
\i j — v'-j. For long lasting contacts, q may not remain parallel 
to the current tangential direction defined by t = v|j/ |v| j | (see, 
e.g,. El) ; we therefore define the corrected = £, ~ n(n • q) 
and introduce the test force 

F t * = -k£-yrm , iJ (2) 

where y, is the coefficient of viscous damping in the tangen¬ 

tial direction (with y, = y„). To ensure that the magnitude of 
the tangential force remains below the Coulomb threshold, we 
constrain the tangential force to be 

P =min^\F n \,\¥ t *\)¥ t */\F‘*\ (3) 

and redefine q if appropriate. 

Cohesive forces are modeled using the approach outlined 
in tB, and are considered to arise from the capillary bridges 
that form when particles get in contact. The functional form 
of this force is given by 

Fb = 27t/?ycos0/(l + 1.05.s + 2.5i 2 ) (4) 

where S = sy/R/V and .? = r (/ - — (di + dj) /2 (taken to be > 0) 
is the particle separation. Here, l/R = 1 /2(1 /d\ + 1 /c/ 2 ) 112011 
(for simplicity we do not account here for polydispersity and 
use d\ = c /2 = 1 in dimensionless units), and V is the volume 
of a capillary bridge between particles. In the present work we 
assume that all capillary bridges are of the same volume. For 
contact angle, 0, we use 0 = 12°, comparable to the value for 
(deionized ultra-filtered) water and (clean) glass 0. For the 
surface tension, y, we use the value corresponding to water, 72 
dyn/cm, scaled appropriately. The critical separating distance, 
s c , at which a bridge breaks is given by 

s c = (1+Q/2)(V 1/3 /R + V 2/3 /R 2 ) (5) 

Here, s c could be thought of as a measure of the strength of 
cohesion; larger s c leads to more pronounced cohesive effects. 

Our simulations are carried out by slowly compressing the 
domain, starting at the packing fraction 0.63 and ending at 
0.90, by the moving walls built of monodisperse particles with 
diameters of size d ave placed initially at equal distances, d aV e, 
from each other. The wall particles move at a uniform (small) 
inward velocity, v c , equal to v ; o = 2.5 • 10~ 5 (in the units of 
dave /T ), or a fraction of it, as we explore the influence of com¬ 
pression speed. Due to compression and uniform inward ve¬ 
locity, the wall particles (that do not interact with each other) 
overlap by a small amount. When the effect of compression 
rate is explored, v c is decreased, or the compression stopped 



(a) F = 0. (b)F = i. 


FIG. 1: (Color online) An example of a reference system for 
different force thresholds at p = 0.9 (see Supplementary 
Material in 0] for animations). 


to allow the system to relax. In order to obtain statistically 
relevant results, we simulate a large number of initial config¬ 
urations (typically 20), and average the results. Due to the 
compression being slow, we do not observe any different be¬ 
havior close to the domain boundaries compared to the rest of 
the domain. 

We integrate Newton’s equations of motion for both the 
translational and rotational degrees of freedom using a 4th or¬ 
der predictor-corrector method with time step At = 0.02. Our 
reference system is defined by N = 2000 polydisperse par¬ 
ticles ( r p = 0.2), with k n = 4 • 10 3 , e„ = 0.5, /j = 0.5, and 
k t = 0.8 k n IB ; the (monodisperse) wall particles have the 
same physical properties. Larger domain simulations are car¬ 
ried out with up to N = 20,000 particles. If not specified oth¬ 
erwise, cohesion is not included. 


HI. RESULTS 
A. Purely repulsive systems 

Figure |TJa) shows an example of the reference system at 
p = 0.90, with the particles color-coded according to the to¬ 
tal normal force, normalized by the average normal force, 
F„/ <F > (we focus only on the normal forces in the present 
work). If the system contains a set of particles in contact that 
connects top/bottom or left/right wall, then there is contact 
percolation. We will also consider force percolation by focus¬ 
ing on the particles sustaining force larger than a given force 
threshold and ask how the percolation properties are influ¬ 
enced by a nonvanishing threshold. As an example. Fig. |TJb) 
shows the same system as in Fig. |TJa) with force threshold 
F = 1. While the system shown in Fig. [TJ a) clearly percolates 
(contact percolation), it is not immediately obvious whether 
the system shown in Fig.QJb) does. 

In describing percolation properties, we use the following 
quantities, all based on averaging over multiple realizations: 
P(p,F), the percolation probability; F p , the percolation force 
threshold, defined by P(p,F p ) = 0.5; and P c (p), the contact 
percolation probability, defined as F c (p) = P(p,0). In addi¬ 
tion, we will use Z, the coordination number, measuring av- 
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(a) The percolation probability, 
P( P,F). 


(b) The percolation force threshold, 
F p , and the coordination number, Z. 


FIG. 2: (Color online) Reference system, averaged over 20 
realizations. 


erage number of contacts per particle; a sharp increase of the 
Z curve is typically associated with the jamming transition, 
see, e.g. 12311 . We note that the listed quantities also depend 
on the number of particles, N, and on the compression speed, 
v e ; this dependence will be discussed later in the paper. For 
the simplicity of notation, we do not include this dependence 
explicitly in the notation. 

Figure |2[a) shows / J (p. F), for the reference system. We 
see that, starting at p « 0.77, there is a percolation transition; 
note that if we vary F and keep p fixed, this transition is rather 
sharp for large p’s and more spread out for p G [0.77,0.81]. 
To describe various transitions that take place as the system is 
compressed, we define: p j, at which jamming, defined here 
as the p at which the Z curve has an inflection point, takes 
place (later in the text we also show that at pj rapid increase 
in pressure (measured at the domain boundaries) occurs, sup¬ 
porting this definition of p /); and p p , at which contact perco¬ 
lation, defined as P c (p„) = 0.5 occurs. Figure [2jb) shows Z 
and F p -, we find from the data shown that pj ss 0.79 (the verti¬ 
cal dashed line in the figure). Note that just below p /, there is 
a strong force network that percolates, as shown by large F p . 
The dominant maximum of F p calls for consideration of an¬ 
other transitional p at which this maximum occurs: however, 
we find that this transition is always sandwiched between p p 
and p j, so we will not discuss it in more details here. 
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FIG. 3: (Color online) Reference system: the percolation 
probability, P c , and Z. 


Figure0a) shows Z and l\ for the reference system. While 


there is some noise in the results, one can still obtain an ac¬ 
curate value for p p ss 0.776. [For this, and all other results 
involving p p and p j, uncertainty of the results is such that the 



(a) fixed compression rate (b) fixed compression speed 


FIG. 4: (Color online) Influence of system size on p p and p j 
for fixed compression rate and for fixed compression speed. 

results are accurate up to three significant digits: for pj we 
use standard error to estimate uncertainty, and for p p we esti¬ 
mate the range over which 0.4 < P, (p ) < 0.6.] Therefore, the 
results for our reference system suggest that p p < p j, and the 
question is whether this finding is robust with respect to the 
changes of the system parameters and of the protocol used. 
Before proceeding, we note that although there are some dif¬ 
ferences between realizations, for all of them we find consis¬ 
tently (for the considered system) that p p and p j differ by a 
nonvanishing amount. 

Regarding the system parameters, we start by discussing 
the influence of polydispersity, measured by r p , and friction 
coefficient, /./. Table [I] shows the results for p p and pj, and 
we observe that both p p and p j are monotonously decreasing 
functions of these two parameters; in particular the results for 
p j are consistent with the ones from literature (see H and the 
references therein). The finding that is perhaps more relevant 
for the present discussion is that the difference between p p 
and p j remains as r p and // are varied. 

Next we discuss the influence of system size; note that this 
issue has been discussed extensively in the context of random 
percolation (see e.g. d). Here, the context is more com¬ 
plicated since the system considered is dynamic, and one has 
to decide on coupling of relevant spatial and temporal scales. 
We have considered two scenarios for the systems of differ¬ 
ent size: one where the rate of the change of p is kept con¬ 
stant, and the one where the compression speed (v c ) is fixed. 
While the details of the results vary depending on the choice 
of the scenario, we find that the difference between p p and 
Pj remains non-zero (and typically increases as a function of 
L) for the both scenarios and for the system sizes defined by 
L = 50, 75, 100, 150: Figure [4] shows the dependence of p p 
and p j behavior on the system size using two aforementioned 
protocols. Figure[4}a) shows results for the fixed compression 
rate; the compression velocity, v c , is increased with L so that 
the rate v c /L is constant. Figure [4jb) shows p p ,pj when we 
keep v c constant as L increases. For both protocols - fixed 
compression rate and speed - we observe increased difference 
between p p and pj as L is increased. 
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0.0 

0.1 

0.2 

0.3 

0.4 

0.5 

Pj 

0.827 

0.812 

0.802 

0.797 

0.796 

0.789 

Pp 

0.815 

0.799 

0.792 

0.784 

0.781 

0.776 



r P 



0.0 

0.1 

0.2 

0.3 

0.4 

Py 

0.804 

0.797 

0.789 

0.7834 

0.782 

Pp 

0.786 

0.784 

0.776 

0.771 

0.766 

v c /vo 


0.0 

0.02 

0.05 

0.1 

1.0 

Pj 

0.798 

0.799 

0.798 

0.792 

0.789 

Pp 

0.798 

0.794 

0.791 

0.786 

0.776 


TABLE I: Influence of p, r p and v c on p p and p j for a 
continuously compressed system (the parameters not 
specified correspond to the reference case). 


Since the reference system is exposed to a nonvanish¬ 
ing compression rate, there is also the question of rate- 
dependence, as already alluded above. To explore this issue, 
we carry out simulations with progressively smaller speed of 
compression, using v c = vo/10, i'o/20 and v'o/50. We find that 
the P c transition becomes sharper as v c decreases, indicating 
that p p is affected by v c ; in general, for a fixed p, the parti¬ 
cles are less likely to percolate for smaller v c and therefore 
p p increases as v c decreases. Both p p and p j are shown in 
Table ID While both p’s increase as v t decreases, the crucial 
finding is that the difference between them becomes smaller 
for slower compression. The question remains whether p p and 
p j collapse to a single value in the limit v c -+ 0. To answer 
this, we consider a modified protocol such that we interject 
relaxation steps in our compression (we reference this proto¬ 
col by v c = 0). More precisely, after compressing the system 
by Sp = 0.001, we check whether there is a percolating clus¬ 
ter. If not, we proceed with compression; if yes, the system 
is relaxed until percolation disappears, and then the system is 
further compressed. We carry out this procedure until such 
p p that percolating cluster does not disappear after relaxation 
(for all considered simulations, the system always percolates 
above p p found using relaxation protocol, or in other words, 
percolation is never found to disappear as a system is further 
compressed). Figure |2b) shows P c and Z for the relaxed sys¬ 
tem, suggesting much smoother and sharper evolution of P c 
through p p . Table IJ shows that for the reference system and 
v c = 0, p p and p j collapse to the same point, within the avail¬ 
able accuracy. We have reached the same finding for the other 
systems listed in Table HI including monodisperse frictionless 
system - while this particular system is known to show differ¬ 
ent behavior due to partial crystallization 0], it still leads to 
p p = Pj. We have also verified that the finding p p = pj still 
holds when different system sizes are considered. 

This finding of collapse of percolation and jamming tran¬ 
sitions appears to be different from the one in flO], where it 
was found that p p and p j differ. The source of the differ¬ 


ence seems to be the use of overdamped dynamics in [0; 
this effect apparently keeps the particles together and leads to 
percolation even for small p’s. We find, however, that, within 



FIG. 5: (Color online) Anisotropy of the stress tensor of 
r p = 0.2, p = 0.5 (squares), r p = 0.0, p = 0.5 (circles), 
r p = 0.2,/j = 0.0 (triangles) and r p = 0.0, p = 0.0 (thick line) 
systems as a function of packing fraction, p. Respective 
jamming transitions, pj , p/, p", p"' , are depicted by a dashed 

line. 


the particle interaction model considered in the present paper, 
based on (constant) coefficient of restitution, e n , the finding 
p p = Pj persists even for very small e n ~ 0, suggesting that 
the finding reported here is robust, within the framework of 
the implemented particle interaction model. 

While the findings obtained in quasi-static limit are of main 
interest, one should note that in the context of particulate mat¬ 
ter, percolation and jamming transitions typically involve dy¬ 
namics, even if very slow one. Close to p j, the relevant time 
scales diverge in the limit of infinite system size, and there¬ 
fore, one could expect that for any sufficiently large system, 
even very slow dynamics may lead to (arbitrarily small) dif¬ 
ferences between p p and p j. Therefore, it should not be sur¬ 
prising if differences are found between p p and p./ for slowly 
evolving spatially extended particulate systems. 

To close our discussion focusing on repulsive systems, we 
discuss whether the implemented compression protocol may 
induce an anisotropy, possibly influencing the results. For this 
purpose, we compute the stress tensor and the distribution 
of the angles of contact between the particles. For brevity, 
we consider here only the compression by i'o- The stress 
anisotropy, x a , is defined by 


<?1 

Ol +09 


( 6 ) 


with Oi ,02 the principal eigenvalues of the Cauchy stress ten¬ 
sor o, specified by o (/ - = 1/(2A)£ c + F,-r,-) as a sum 

over all inter-particle contacts Ck for all particles p; (wall par¬ 
ticles as well as the contacts of interior particles with the wall 
particles are not included here). Here, A is the total area of the 
system, rp rj are the x and y components of the vector point¬ 
ing from the center of particle p towards the particle contact 
q. +,, Fj denote the x, y components of the interparticle force 
at the contact c>. 

Figure [5] shows x a as a function of p. We depict jamming 
transitions, p j, p'j,p" and p'J' by dashed lines for p = 0.5, r p = 
0.2 (reference system), p = 0.5, r p = 0.0, p = 0.0, r p = 0.2 and 
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cn = 2 — 6 contacts (circle, square, 
delta, gradient and thick line, 
consecutively). 



FIG. 6: (Color online) Distribution of the angles of contacts 
for the reference (a), (b) and r p = 0.0, ,u = 0.0 system (c), 
(d). In these polar plots, the azimuthal coordinate, (j). 
corresponds to the angle between the line connecting the 
centers of contacting particles and the +x axis, and the radial 
one to the probability of observing given (j). 


ju = 0.0, r p = 0.0, respectively. While far below the jamming 
(and percolation) transitions, the anisotropy measured by x a 
may be present, close to p p and pj, x„ C 1 for all systems 
considered, showing that the systems are essentially isotropic 
for the packing fractions of relevance here. Above jamming 
points, T a is even smaller. 

Figure [6] shows the distribution of contact angle, (j), for 
the reference system, the parts (a), (b), and for the ju = 0.0, 
r p = 0.0 system, the parts (c), (d). Most importantly, this fig¬ 
ure shows symmetric distribution of 4>’s. In addition, by com¬ 
paring the results of the reference case with the ones obtained 
for monodisperse frictionless, we also observe the influence of 
partial crystallization on the latter, for large packing fractions. 


B. Cohesive systems 

Here, we discuss the effect of cohesion on percolation and 
jamming. We have considered few different ‘strengths’ of 
cohesion (specified by the distance, s c ), at which capillary 
bridges break; for brevity here we present results only for 
‘weak’ cohesion, specified by small distance at which cap¬ 
illary bridges break, s c « 0.0028 -C 1 (see SecITili. We focus 
on the relaxed reference system. Figure [7] a) shows that the 
percolation transition occurs very close to (the starting value) 
p = 0.63. The Z curve remains at high values for all con¬ 


FIG. 7: (Color online) (a), (b) Cohesive relaxed system, (c), 
(d) Pressure on the walls, II, and Z. Dashed lines correspond 
to p/ (and to p p in (c)); in (a), p p is shown by dotted line. 


sidered p’s, but we note that there is a kink in the Z curve 
at p ss 0.783. The kink and consecutive increase of Z sug¬ 
gest that the system undergoes a transition. To verify that 
this transition corresponds to p j, we consider the pressure on 
the system walls, II. Figure |7Jc and d) shows this pressure 
(force/length, in dimensionless units) for both the reference 
system, and for the cohesive one. We see that for the refer¬ 
ence system an increase of II occurs at pj (inflection point of 
the Z curve). FigurcfTJd) shows that an increase in II and the 
kink in the Z curve occur at the same p = pj = 0.783. 

Clearly, the difference between pj and p p is significant for 
the considered cohesive system, consistently with the earlier 
work III. As expected, we find similar results for the systems 
characterized by larger s c (results not shown for brevity). The 
strong influence of weak cohesion on the p p and p j suggests 
that for any non-vanishing cohesion, one would find differ¬ 
ences between p/j and p j, with this differences disappearing 
only in the limit of s c —>• 0. As soon as there is no attractive 
force, the difference between p p and p j vanishes even in the 
limit of inelastic collisions, e„ —> 0. 

One may ask about the origin of the ‘kink’ in the Z curves 
for the cohesive system. An intuitive explanation is as fol¬ 
lows: as compression starts, the particles immediately get in 
contact, form mini-clusters (consisting of a small number of 
particles), leading to rather large Z; due to the presence of co¬ 
hesive forces, relaxation does not lead to breakup of the exist¬ 
ing contacts. Therefore, as long as p is small, the mini-clusters 
do not break; as p grows, however, collisions start separating 
particles, leading to breakup of the mini-clusters and decreas- 
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mg Z. At some point, when p becomes sufficiently large so 
that all particles are effectively in contact, Z starts growing 
again, and at the same p, II starts increasing. To support this 
description. Fig. [7jb) shows the number of particles ( N ) with 
2,...,6 contacts ( cn ). We observe that as p/ is approached 
from below, the cn = 4, 5 curves have negative slope, sug¬ 
gesting breakup of the clusters (this breakup is presumably 
also partially responsible for the positive slope of c n = 2, 3 
curves for the same values of p); at py these trends reverse. 

IV. SUMMARY AND CONCLUSIONS 

Percolation and jamming transitions of evolving particulate 
systems are non-trivial. We find that these transitions for re¬ 
pulsive particles interacting by a commonly used interaction 
model coincide for quasi-static systems; this finding, together 
with the results reported in fToll , where these transitions are 
found to differ for particles following overdamped dynam¬ 


ics, suggests that the considered transitions may be influenced 
significantly by the type of interaction between the particles. 
Furthermore, our finding is that any, even very slow dynamics 
may lead to the differences of the packing fractions at which 
percolation and jamming occur. Therefore, in particular close 
to jamming, a careful exploration will be needed in order to 
distinguish the effects due to dynamics and due to, e.g., the 
type of interaction between the particles. In the same vein, we 
are also finding that even minor cohesive effects have a strong 
influence in particular on percolation transition. 

We hope that the present results will encourage carrying out 
careful experiments that will quantify further the predictions 
regarding the influence that dynamics, cohesion, and the na¬ 
ture of particle interaction have on percolation and jamming. 
Our own research will continue in the direction of exploring 
the effects of jamming and percolation in three spatial dimen¬ 
sions. 
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